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The feasibility of path integral Monte Carlo ground state calculations with very few beads using a 
high-order short-time Green's function expansion is discussed. An explicit expression of the evolution 
operator which provides dramatic enhancements in the quality of ground-state wave-functions is 
examined. The efficiency of the method makes possible to remove the trial wave function and thus 
obtain completely model-independent results still with a very small number of beads. If a single 
iteration of the method is used to improve a given model wave function, the result is invariably a 
shadow-type wave function, whose precise content is provided by the high-order algorithm employed. 

PACS numbers: 31.15.Kb,02.70.Ss 
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Quantum Monte Carlo path integral calculations pro- 
vide a powerful approach to many-body physics, both 
at zero QUI and finite temperature [4(. They rely 
on using the classical action of the system in imaginary 
time to define a path distribution probability, following 
Feynman's approach. For systems with Bose statistics, 
the probability is positive and the only approximation 
involved in the calculations comes from the discretiza- 
tion of time, which prevents an exact evaluation of the 
action. In path integral calculations, chains with large 
numbers of time slices or beads are just a realization of 
these ideas [j| . If an accurate evaluation of the action for 
paths with large time steps were possible, the complexity 
of path integral calculations, and in particular ergodicity 
issues, would be greatly reduced. 

For ground-state calculations, the path integral ap- 
proach (PIGS) [3, S i provides a method to systemati- 
cally improve a trial wave function, by repeated applica- 
tion of the evolution operator in imaginary time, which 
eventually drives the system into the ground state [f|. 
Here again, if an accurate implementation of the evolu- 
tion operator for large time steps were available, it could 
be used to build variational wave functions of very high 
quality. 

In this letter, we focus on high-order short-time expan- 
sions of the Green's function, which have been the subject 
of a series of works in recent years @, 8, 9]. This progress 
has led to a deep understanding of their properties and 
provided various extremely accurate decompositions of 
the evolution operator U — exp(e H), with H = T + V , 
expressed as products of unitary operators of the basic 
components U — exp(e T) and U — exp(e V). An im- 
portant result is that the use of a finite time step e can 
be renormalized in the potential. This means that for 
a range of time steps e, one is in practice virtually free 
from finite time step errors. 

We make use of the decomposition schemes proposed 
by Chin in order to evaluate the feasibility of PIGS 

calculations with high-order propagators using very few 
beads (Nf,). To this end, we revisit the ground state 



of bulk superfluid 4 He. We also explore the possibility 
of performing ground-state calculations of boson systems 
without any model wave function, i.e., to take as a start- 
ing point \P = 1 and rely in the propagator's quality 
to build the actual wave function exclusively from the 
Hamiltonian. It is also interesting to note that if one 
does make use of a trial wave function, the application 
of Chin's evolution operator [ij [1(3] produces a much en- 
hanced model that actually incarnates a shadow wave 
function [6( . In this way, Chin's analysis provides a deep 
understanding of the success of shadow-like wave func- 
tions [ill, Ell and sheds light onto the actual mechanisms 
leading to this kind of wave-functions. 

Decompositions of the evolution operator preserving 
unitarity in the form 

N 

exp(e (T + V)) = TJ expfeef) exp( Vi eV) + 0{e n+1 ) , 

i=l 

(1) 

arc the starting point which have culminated in the al- 
gorithms due to Chin that we test here. By means of 
a proper selection of the factorization coefficients {U} 
and {fi}, any given order can be achieved (the result- 
ing expression is nth order since the effective Hamilto- 
nian is then T + V + 0(e n )). However, if some of the U 
factorization coefficients are negative, involving diffusion 
processes backwards in time, the algorithm cannot be 
used in the context of quantum Monte Carlo calculations. 
Fourth order algorithms, for which Chin has worked out 
a complete characterization happen to be somehow 
unique, as for only forward decompositions, the highest 
order that can be achieved is four. These developments 
result from a careful use of the expansion 

N 

Y[ expfaef) exp(w l eV r ) = exp(ee T T + ee v V + e 2 e TV (2) 

8=1 

x [f , V] + e 3 e TTV [f, [f, V}} + s 3 e VT v[V, [f, V}]) + ... 
which allows to keep only the simplest term [V, [T, V]], 
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which produces a sort of renormalized potential 

h 2 N 
i=i 

with Fi = Yljjii^iV{ r ij)- A. remarkable characteristic 
of the resulting expansions is that they lead to a con- 
tinuous family of 4th order algorithms characterized by 
one free parameter. It is interesting to note that with a 
proper choice of that parameter {t\ in Eq. [4]) one is able 
to fine tune the leading error term of the propagator to 
any desired value, including changing sign in a controlled 
way. Thus, one is enabled to minimize the value of the 
4th order error coefficient, and in doing so, even try to 
cancel, to the largest possible extent, the contributions 
of the next orders. In practice, this means that a partic- 
ular value of the free parameter entering in the algorithm 
produces exceedingly accurate and stable results. A dif- 
ferent strategy to improve the order of the propagator 
is the use of a multiproduct expansion with a controlled 
violation of positivity [HI]. 

Amongst the various decompositions proposed in Ref. 
Q, the particular scheme chosen in this work is 

e -TH _ e -VQT V g-tlT f g-VlT W g-*2T f g -«lT W g-tlT f g-l>OT V 

(4) 

with W = V + (uq/v^t 2 [V, [f , V]}. The various factor- 
ization coefficients are all dependent on the single free 
parameter t\ The range of possible values for tiis 
< t\ < 0.5 and experience shows that the optimal value 
is nearly independent of r = e/Nb- A similar decomposi- 
tion has been recently used in path integral Monte Carlo 
(PIMC) simulations at finite temperature showing high 
accuracy 

The operator splitting (2J allows for an estimation of 
the Green's function G(R, S,r) = (R\e- rH \S), which in 
turn provides the mechanism for building accurate wave 
functions starting from a trial wave function, ^f(R) = 
J G(R, S, r)^ m (S) dS. More explicitly, 

*(R) = J e -»o^(fl) e -^ e - V w(S l)e -%# (5) 
x e - VlT w(s 2 ) e - { % t s ^ e -v rv(s 3 ) y m (s 3 ) dSidS 2 dS 3 . 

Actually, Eq. sets the grounds for building wave 
functions consistent with the decomposition of the prop- 
agator ((4}. It has an intuitive content too: it states that 
the actual value of *B(R) should be taken as a weighted 
average of neighbor shadow values 'I'(S'), with a weight 
given by a precise combination of exponentials of the 
(finitc-timc-step renormalized) potential. Moreover, it 
is also very suggesting to consider the double possibility 
that Eq. [5] offers, either as a powerful enhancement of 
an a priori model wave function ^ m (i?) or directly as 
a tool to generate the wave function directly from the 
Hamiltonian. 



Equation \E\ can also be written in terms of relative 
rather than absolute auxiliary coordinates, 

e -v rV{R) e - 4Dt ' 1T e -firVy(-R+Si) e - 4J>f 2 2T 
x e ~v 1 rW(B.+S 1 +S 2 ) e - {R+S lt,^~ S3)2 e -v rV(S 3 ) 

x * m (i? + Sx + S 2 + S 3 ) dS!dS 2 dS 3 (6) 

Whilst Eqs. (J5]) and (HJ) are equivalent, they lead to dif- 
ferent estimators for the kinetic energy. The action of the 
kinetic operator on ^(i?) in Eq. © involves derivatives 
of \I/ m , in contrast to what happens when the prescription 
given in Eq. ([5]) is used. 

It is possible to use the propagator G(R, S, r) in order 
to improve the quality of any given trial wave function, 
with reduced variance in direct proportion to its quality. 
In fact, Eq. ([6|) satisfies the principle of zero variance: as 
the trial wave function ^ m approaches the exact ground- 
state wave function, the propagation time r can be con- 
tinuously tuned to smaller values, with the S distribution 
approaching a Dirac delta and ^(R) and its derivatives 
approaching the exact ones. 

In order to test the accuracy of the method (|5l6p . 
we have applied it to liquid 4 He, a classical benchmark 
in quantum many-body physics. The calculations have 
been carried out at the experimental equilibrium density 
p = 0.365 a~ 3 (a = 2.556 A) with TV = 64 atoms and us- 
ing the HFD-B(HE) Aziz potential [HI, which has proven 
to be highly accurate in the description of the experimen- 
tal equation of state As a trial wave function, wc 
use a simple Jastrow form based on the McMillan model, 
i£ m (R) = Y[ %] exp[-0.5(6/r y ) 5 ] with b = 1.20 a, and two 
additional variational parameters, t\ and r ©. After a 
quick search we found the optimal value t\ = 0.35 which 
was kept fixed for the rest of the calculations. The de- 
pendence of the variational energy on the remaining vari- 
ational parameter r is shown in Fig. [TJ where the varia- 
tional character of the calculation is clear. The horizontal 
axis corresponds to the r parameter in K _1 units and the 
vertical axis represents the total energy per particle in K. 
The data displayed as empty circles is the variational en- 
ergy for the wave function based on the high-order action 
([5])- The position of the minimum is a compromise be- 
tween a large r value suitable for a large suppression of 
the excited components present in \l/ m (i?), and a small 
one suitable for a proper behavior of the variational wave 
function (JSJ), which relics itself on a short-time expansion. 
The minimum is located at r = 0.025 and the variational 
energy obtained is E = —7.10 K, only ~ 0.2 K higher 
than the exact value. 

Equally impressive is the data displayed as empty dia- 
monds, which is a variational calculation using the same 
propagator (O acting on ^> m (R) = 1, i.e. it is a varia- 
tional calculation in which only the Hamiltonian and the 
statistics are used. We see in this case that the compro- 
mise between a large suppression of the non-ground-state 
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e (K- 1 ) 


N b 


McMillan 
E/N (K) V/N (K) 


Semi-classical 
E/N (K) V/N (K) 


E/N (K) 


1 

V/N (K) 


0.01 


1 


-6.860(58) 


-20.924(30) 


-3.878(76) 


-17.246(44) 


-2.81(9) 


-17.009(52) 


0.02 


2 


-7.175(51) 


-21.106(37) 


-6.234(76) 


-20.285(44) 


-5.73(12) 


-20.148(42) 


0.04 


4 


-7.268(44) 


-21.413(35) 


-7.121(67) 


-21.284(43) 


-6.93(11) 


-21.186(39) 
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-7.303(35) 


-21.538(44) 


-7.306(64) 


-21.547(40) 
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0.08 
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-7.299(41) 


-21.534(38) 


-7.290(55) 


-21.583(39) 
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-7.25(8) 


-21.381(20) 


0.10 


8 










-7.33(8) 


-21.498(12) 


0.12 


10 










-7.34(9) 


-21.524(10) 



TABLE I: Total and potential energies per particle as a function of the imaginary time e and number of beads Nt using different 
models for the initial wave function ^ m (i?). 



components and a small r parameter has been shifted to- 
wards a higher value r = 0.04, as could be expected. 
A third model for ^> m (R) consists in a semi-classical 
approximation, ^ m _ sc (R) = JJ i<3 - exp[-EV SI (rij], with 
V sr (r) — 1/r 12 corresponding to the r-dependence of 
the Lennard- Jones potential around the core. This third 
model does not contain any free variational parameter 
since e is the total imaginary time of the propagator. 
The results obtained with 'fm-sclV) (open squares) only 
improve slightly the variational energy obtained with 
V m {R) = l. 

It is also possible to apply several times the propagator 
to ^> m (R) in order to obtain better variational results 
and finally projecting out any excitation present in the 
model wave functions. This is shown in Fig. [T] with a 
filled circle point, which corresponds to the propagator 
applied Nb — 4 times onto the Jastrow-McMillan wave 
function ^f m (R) for a total time propagation e = 0.04 
K _1 . The same figure shows with a filled diamond the 
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FIG. 1: Variational energies obtained with a single time step 
r and the proposed wave function. Open circles, squares, and 
diamonds correspond to *& m (-R) of McMillan type, semiclassi- 
cal model and 1, respectively. The corresponding filled sym- 
bols stand for the asymptotic values using the three models 
and more than one time step. 



result of acting with Nb = 10 propagators on ^ m (R) = 1 
for a total time propagation e = 0.12 K _1 . Between these 
two points, the semiclassical model achieves convergence 
very fast with Nb — 4 and time e = 0.06 K _1 . 

The high-order Green's function can also be used to 
obtain the total energy following the standard procedure 
of evaluating the local energy of the trial wave function 
^f m (R) at one end of the chain, and its high accuracy 
translates again in the need for a very small number of 
beads until convergence is achieved. The results obtained 
for the total and potential energies per particle are re- 
ported in Table U as a function of the total imaginary 
time s. The potential energies are calculated in the cen- 
ter of the chain to remove any possible bias coming from 
^ m (-R); the total energies are estimated in the extreme, 
except for the case <3/ m (-R) = 1 where they are sampled 
in the center. The number of beads for a given time e 
is determined simply by the requirement that doubling 
its number (and simultaneously halving the propagation 
time per bead r) does not have any effect, which turns 
to be equivalent to keep the propagation time per bead 
below t = 0.15 K _1 . One would expect that this regime 
corresponds to keeping the finite time step error just be- 
low the detectable level. This is in accordance with Fig. 
[TJ where we see that a propagation time per bead near 
0.15 K _1 is close to the maximum value before the time 
step error starts to become apparent by bending upwards 
the variational energy curve. When both requirements: 
t = e/Nb to be small enough (t < 0.15 K _1 in our case) 
and Nb large enough are met, the energy becomes inde- 
pendent of both e and JV& and a good estimation of the 
ground-state energy is achieved. 

The results of Table U show that the convergence 
is quickly achieved with only a few number of beads: 
A5 = 4 for the McMillan and semiclassical Jastrow fac- 
tors, and Nb = 8 for ^ m (R) — 1. Concerning the con- 
vergence for the potential energy, one can see in Table U 
that the exact (asymptotic in e) value is independent of 
the trial wave function and its value is reached with only 
Nb = 4. It should be taken into account however that, as 
in any PIGS method, the total length of the chain corre- 
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FIG. 2: Two-body radial distribution function g(r). The solid 
and dotted lines correspond to present results using for ty m {R) 
a Jastrow-McMillan factor or ^ m (R) — 1, respectively. The 
dashed line is the experimental data from Ref. [TtJ 
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has been developed for path integral Monte Carlo (finite 
temperature) and that we have extended to PIGS. One 
of the main advantages of this method is that p(r) comes 
properly normalized, and thus eliminates any uncertainty 
introduced by the a posteriori normalization factor. In 
Fig. [3J results for p(r) obtained using different trial wave 
functions ^ m are shown. As one can see, the results are 
statistically indistinguishable and predict a condensate 
fraction no — 0.080(2), in nice agreement with a recent 
PIMC estimation at T = 1 K (n = 0.081(2)) 0. 

Summarizing, the high-order short time expansion of 
the Green's function presented in this work enables the 
possibility of performing variational calculations of very 
high-quality on systems for which no model of wave func- 
tion is known. This is illustrated in the case of liquid 
He at equilibrium density, where the propagator pro- 
vides for Nb = 1 and *ff m (R) = 1 a variational energy 
E/N — —6.20 K, and converges to the exact value al- 
ready with Nf, — 8. When used to improve a Jastrow- 
McMillan wave function, it produces a shadow-like vari- 
ational wave function for Nb = 1 with a variational 
ground-state energy E/N = —7.10 K, only ~ 0.2 K 
higher than the exact value. Repeated application of the 
propagator leads to variational results which are asymp- 
totically exact for values as low as Nb = 4. The prospects 
for future work are promising, since this opens the road 
to being able to obtain results for systems whose ground- 
state wave function is poorly known or even unknown. 

We wish to thank stimulating discussions with Siu 
Chin. Partial financial support from DGI (Spain) Grant 
No. FIS2008-04403 and Generalitat de Catalunya Grant 
No. 2005SGR-00779 is also akcnowledged. 



FIG. 3: One-body density matrix p(r). Filled circles and open 
squares stand for PIGS results using for ^ m {R) a Jastrow- 
McMillan factor (N b = 5) or * m (i?) = 1 (N b = 10), respec- 
tively. The error bars are smaller than the size of the symbols. 



sponding to |W m (i?)| is twice that of ty m (R), and that 
one r propagator ([5]) involves three internal shadows. 

Unbised (pure) estimations of operators O other than 
the Hamiltonian can only be calculated in the center of 
the chain, (6) = Af- 1 (^ m \G(s/2)OG(e/2)\^ m ). This 
holds, for instance, for the potential energy reported in 
Table U and the two-body radial distribution function 
g(r) shown in Fig. O The present PIGS results for gjr) 



show an excellent agreement with experimental data 17 1 
for the trial wave functions used in this work. Even when 
^m{R) — 1 the result is the same, the calculation requir- 
ing only a few more beads. 

Another relevant function that can be computed in an 
unbiased way is the one-body density matrix p(r), whose 
asymptotic limit is the condensate fraction. The cal- 
culation of p{r) has been carried out by incorporating 
worm movements in the sampling [l8|], a technique that 
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